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Abstract 

A simple and efficient variational method is introduced to accelerate the conver- 
gence of the eigenenergy computations for a Hamiltonian H with singular potentials. 
Closed-form analytic expressions in TV dimensions are obtained for the matrix ele- 
ments of H with respect to the eigenfunctions of a soluble singular problem with two 
free parameters A and B . The matrix eigenvalues are then optimized with respect 
to A and B for a given N . Applications, convergence rates, and comparisons with 
earlier work are discussed in detail. 
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1. Introduction 

Sources of ongoing interest [1-66] in the study of singular potentials are at least three- 
fold: (i) regular Rayleigh-Schrodinger perturbation theory can fail badly for such 
potentials, and (ii) in physics, it is quite common to encounter phenomenological 
potentials that are strongly singular at the origin, and (hi) because of the intrinsi- 
cally interesting mathematical problems that arise in their study. A specific family 
of singular quantum Hamiltonians that has found widespread application in many 
areas of atomic, molecular, nuclear physics, are the so-called anharmonic singular 
Hamiltonians given by 

where r e (0, oo) and a n and A„ are positive real numbers. A particular and important 
special subclass of (1.1) that has been a subject of intensive studies is the set of spiked 
harmonic oscillator Hamiltonians [9-3] 

H = -^+r 2 + ± a>0 ,A>0 (1.2) 

acting in the Hilbert space L 2 (0,oo) with eigenf unctions ip e L 2 (0, oo) of H satisfying 
the Schrodinger equation -ip" + (r 2 + \r^ a )ip = Eip with ^(0) = 0. The function V is 
an eigenfunction corresponding to the eigenvalue E and the condition ip(Q) = is the 
Dirichlet boundary condition we impose. The coupling parameter A determines the 
strength of the perturbative potential and the positive constant a represents the type 
of singularity at the origin. Thanks to the pioneering work of Detwiler and Klauder 
[8] and of Harrell [9] , remarkable progress has been made over the past three decades 
in the field of spectral calculations for the spiked harmonic oscillator Hamiltonian 
(1.2). Other interesting subclasses of (1.1) which have been used in atomic, molecular 
and nuclear physics are the anharmonic singular Hamiltonians [47-61] 

d 2 

H = -— + ar 2 + br~ 4 + cr" 6 , a>0, 6>0, c> (1.3) 
dr z 

and the positive-parameter singular even-power Hamiltonians [62-66] 

rj d 2 2 a 2 a 3 a 4 , . 

Several numerical and analytical (both variational and perturbative) techniques are 
available in the literature for the exact and approximate eigenvalue calculations for 
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these families of Hamiltonians. Some of these techniques were devoted to particular 
classes of singular Hamiltonians and some others were restricted to specific values of 
the potential parameters. It is of great interest to have a successful method valid to 
study all of these classes without the need for major adjustment as one goes from 
class to class. The purpose of this article is to provide such technique and to present 
a rigorous variational approach for the accurate calculation of the energy levels of the 
singular Hamiltonians (1.1). Our method is very simple and yet accurate enough to 
determine the entire spectrum in arbitrary dimensions of anharmonic singular Hamil- 
tonians (1.1). It can be viewed as an extension of the earlier variational approach to 
the study of the spiked Harmonic oscillator potentials [44]. The principal ingrediants 
are (a) the use of a special basis comprising the exact solutions of a singular problem, 
(b) the determination of exact analytical expressions for the matrix elements with 
respect to this basis, and (c) the retention of two parameters from the basis problem 
that can be used as additional variational parameters. These optimizations are made 
possible because of the progress in the establishment and simplification of closed-form 
analytic expressions for the matrix elements. 

The article is organized as follows. In Section 2, we introduce our variational 
technique: the method is discussed in a general setting with no reference to particular 
application. In Section 3, we introduce the Gol'dman and Krivchenkov Hamiltonian 
as a solvable model. Thereafter, we use its eigensolutions to compute the matrix 
elements for more general singular operators r~ a , a > . Closed analytical expressions 
in terms of single finite sums are obtained for the matrix elements of the power- 
law potentials r q ,q = 2,4,6, In Section 4, the applications to singular potentials 

are discussed. Special attention is paid to spiked harmonic oscillator (1.2) where we 
compared our results with other techniques available in the literature. Thereafter, we 
study various higher-order anharmonic singular potentials, such as (1.3) and (1.4), 
using the approach discussed in Section 2. The convergence problem of the variational 
approach is studied in some detail. Comparisons with different methods for special 
classes of anharmonic singular potentials are also studied. 

2. Methodology 

In this section, we develop a detailed variational method for studying the family 
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of Hamiltonians (1.1). The variational function is taken as a linear combination of 
orthonormal functions of an exactly solvable model which itself has a singular poten- 
tial. Let {E n (Q,),ip%} be the eigenvalues and the eigenfunctions of an exactly solvable 
Hamiltonian H n = -A + v n (r) acting on separable Hilbert space L 2 {D,d^) such that 



n ' 



I! V>£ \\l= <V£lV£> = / mr)?dKr) = 1. 

Here u a (r) is taken as a function of r and depends on the parameters in the set £1 . 
Let H A = -A + V A (r) be the quantum Hamiltonian under investigation, where A is a 
fixed set of parameters, and we let e(A) be the associated exact eigenvalues of H A . 
By writing the Hamiltonian H A in the extended form 

H A = H n + V A (r)-v n (r), (2.2) 

we have that if the integrals 

<^\V A (-)\^> = J i%V A ^ are finite, (2.3) 

D 

then the matrix elements of the Hamiltonian H A can be written, for m,n = 0, 1,2, .. . 
as 

H A n ee <1%\H\tf> = E mn (Q)5 mn + <til\V A (-)\^> - <^\v%)\^>, (2.4) 
or in compactified form as 

hL = <€\ A|^> + <^\v A (-)\^>. (2-5) 

Our variational technique is based on forming trial wave functions from a linear 
combination of D -orthonormal functions Vn( r ) , n = 0, 1,2, 1 

D-l 

*(r) = Yl c ^n(T) (2.6) 

n=0 

The linear parameters c„ that optimize the energy are determined by the following 
system of equations 

15-1 

J2(H A n -e(A)5 mn )c n = 0, ro = 0,l,2,...,£>-l, (2.7) 

n=0 
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The necessary and sufficient condition for a nontrivial solution of (2.7) is the vanishing 
of the secular determinant 

det|fl£„ - e(A)<y mn | = 0, (2.8) 

The condition (2.8) yields upper bounds to the exact eigenvalues e(A) by means of 
the inequality 





/ H oo 


H 01 


H 0D-1 \ 




e(A) < min diag | 


rrA 
-"10 


H& . 


H 1D-1 


(2-9) 






H D-11 ■ 


fM 1 





where H A n = <^\H A \^> = H% n + v£ n (r) - v% n (r) , for fixed A, are functions of 
the parameter set il . Note that the equality in (2.9) holds if V A (r) = v n (r) ^-a.e. 
in r . The computation of the right hand side of (2.9) requires a diagonalization of 
the matrix over the D -dimensional subspace spanned by orthonormal functions 
followed subsequently by a minimization over the parameters fl . The advantages 
of this method are: (i) only a few matrix elements are needed to achieve accurate 
bounds to the eigenvalues; (ii) The exact eigenvalues are approached monotonically 
as D is increased; (iii) the minimization over a set of parameters accelerates the 
convergence of the energy bounds more rapidly than any standard minimization over 
a single variable; (iv) the optimization of (2.9) for DxD matrix gives upper bounds to 
the energy eigenvalues of the lowest D -states; (v) the diagonalization of such a matrix 
also produces the coefficients required for the corresponding eigenvectors determined 
variationally. 

3. Exactly solvable model and associated matrix elements 

The accuracy and the computational simplicity of the variational method depends 
greatly on the analytic structure of the wave functions that we use, in particular, 
their behaviour in the neighborhood of the singularity. Many different forms of trial 
wave function have been explored in the literature to solve the spiked harmonic oscil- 
lator problem (1.2). The rate of convergence for a variational calculation depends on 
the ability of the basis functions used in the variational calculation to approximate 
the behavior of the exact wave function in the neighborhood of the singularity. Re- 
cently, Hall et al have pointed out the advantages of basing the variational analysis 
of singular potentials on an exact soluble model which itself has a singular potential 
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term. They have suggested and used trial wave functions constructed by means of the 
superposition of the orthonormal functions of the exact solutions of the Gol'dman 
and Krivchenkov Hamiltonian 

*o = -£s+Br* + ±. (3-D 

It was shown that this orthonormal basis serves as an effective starting point for the 
variational analysis of the Hamiltonian (1.2). In this paper we use these solutions of 
H to provide systematic variational solutions for the singular Hamiltonians (1.1). 
The Gol'dman and Krivchenkov Hamiltonian (3.1) is one of the few tht admit exact 
analytical solutions. The Hamiltonian is the generalization of the familiar harmonic 
oscillator in 3-dimension -d 2 /dr 2 + Br 2 + 1(1 + l)/r 2 where the generalization lies in 
the parameter A ranging over [o, oo) instead of only values determined by the an- 
gular momentum quantum numbers I = 0,1,2, The background on Gol'dman and 

Krivchenkov potential V(r) = Br 2 + Ar~ 2 relevant to the following discussion can be 
found in [44]. In particular, the energy spectrum of the Schrodinger Hamiltonian H 
is given, in terms of parameters A and B , by 

£„ = 2/3(2n + 7 ), n = 0,l,2,..., (3.2) 



in which (3 = \[~B and 7 = 1 + J A + 4 and the normalized wavefunctions are 



Mr) = (-irJ^^fr^ie-^^i-n^Pr 2 ). (3.3) 



Here 1F1 is the confluent hypergeometric function 

iFi(-n; b; z) = y —— — - , ( n -degree polynomial in z) 
k =o i - b > kk - 

and the shifted factorial (a) n is defined by 



(3.4) 



(o)o = l, {a)„ = a(a + l)(o + 2) . .. (a + n- 1), for n = 1, 2, 3, . . . . (3.5) 

may be expressed in terms of the Gamma function by (a)„ = T(a + n)/T(a), when a is 
not a negative integer -m , and, in these exceptional cases, (-m)„ = if n > m and 
otherwise (-m)„ = (-l) n m!/(m - n)\. 

An important observation regarding the solvable model H is the existence of the 
A -term which has the dimensions of kinetic energy such as the term that appears 
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in higher-dimensional systems. This observation allow us to extend (3.2) and (3.3) 
to the exact solutions of N— dimensional Gol'dman and Krivchenkov Hamiltonian, 
namely, 

-£ + 4*±£^ + fl * M>O.B>0> (3.6, 

where A = \{M — 3) , and M = N + 21 . The exact solution of (3.6) can be easily found 
by replacing A in (3.2) and (3.3) with 

A-»A(A + 1) + A (3.7) 

This particular observation can be extended to any 3 -dimensional exact solvable 
quantum model. Indeed, if Schrodinger's equation can be solved for arbitrary angular 
momentum number / , then the extension to the N -dimensional case can be obtain by 
replacing I with A . It should be also noted that N and I enter into the Hamiltonian 
(3.6) in the form of combination N+21 . Hence, the energy for a spherically-symmetric 
potential V(r) are the same as long as M is not altered. We now summarize the 
exact eigenvalues of N -dimensional Schrodinger equation with the Gol'dman and 
Krivchenkov potential as 

E% = 2p(2n + lN ), n,l = 0,1,2,..., (3.8) 

where (3 = \f~B and -f N = 1 + ^ A + (A + \) 2 , while the exact eigenf unctions are given 
explicitly by 



= 2 -^^r^-h-^\ Fl ( ^ |/?r 2 ) , (n, I = 0, 1, 2, . . . , N > 1). (3.9) 

In the next section all our results are forumlated in arbitrary dimension N > 1. 
3.1 Matrix elements of the singular operator r Q : 

The effectiveness of the variational method relies on finding a basis that allows for easy 
calculation of the matrix elements of the given Hamiltonian. An important advantage 
of the orthonormal wavefunctions (3.3) is the existence of closed-form formulas for the 
singular potential integrals <i) m \r~ a \ip n > ■ These closed form expressions are achieved 
by means of the following identity: For m and n non-negative integers and 27 > a 



00 

/ ^-"-'e-<*>. ( 7 k) ( 7 1-* 2 )* - (""■- 





(f) n r( 7 -f) „ ^-m, 7 -f,l-f 



7,1- f -n 



1 

(3.10) 
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where the Clausen hypergeometric function 3 F 2 is defined by the series representation 

( —m,a,b A (— m)k(a)k(b)k , j i • i\ 

'( c,d 1 )=^ o (cUd) k kl ' (-"degree polynomial). 



3-^2 



The proof of this identity and some relevant integrals can be found in [44]. Thus the 
matrix elements r m % = <ip n \r~ a \ip n > of the singular operator r~ a have the explicit 
forms 



= ( 1 y t +m ;? # (f)" r ('T-f) ,/ (7)n(7)m ^ ( ~m^f - § ,1 - § 



1 . (3.11) 



(7)« T( 7 ) V n!m! 7, 1 - f - « 

In the case of a being a non-negative even number (a = 2,4,6,...), the Clausen 
hypergeometric function 3 F 2 in (3.11) may be looked upon as a polynomial of degree 
f - 1 instead of an m -degree polynomial. This is, of course, not the case for < a ^ 
2,4,6,... in which case the numerical computational would have to be done directly 
using the expression (3.11). For n > m and a = 2,4, 6,... we have by means of the 
series representation of the hypergeometric function 3 F 2 that 



3 r 2 



-(f-l),7 
7,1-f 



2 ' 

- n 



- rn 



f- 1 

E 

s=0 



(_m) s ( 7 -§),(!-§) 



2 ^ s 



*!(7).(l-f -n) 4 



(3.12) 



As a result, the matrix elements (3.11) further simplify into the closed form expres- 
sions immediately appearing. These are most suitable for computational purposes as 
for the case of 7 > 1 and a = 2 , we indeed have 



I 1 \m+n P l n\(i) m 
\ L ) 7-1 V rn\(~t)„ 

_§_ 

7-1 



(-1) 



m+n I m\(j 



7— 1 y 711(7)^ 



if n > m, 
if n — m, 
if n < m. 



(3.13) 



On the other hand, for 7 > 2 and a = 4, we have from (3.11) and (3.12) that 



(-ir+",3 2 



'(7), 



7( 7 -i)(7-2) V [7(n - m + !) + 2m l if n>m, 

7(7-^7-2) ^ + 

(-l) m+ "/3 2 
I 7(7-l)(7-2) 



(3.14) 



y™^[ 7 ( m -n+l) + 2n] if n<m. 
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We also point out, for 7 > 3 and a = 6, Eq.(3.12) lets us deduce 

(-l) m + "/3 3 /"!( 7 )r, 



2( 7 +l) 7 (7-l)(7-2)(7-3) V "»!( 7 )„ ~ 

[(2 + n)(l + 71)7(7 + 1) - 2m(l + n)( 7 - 3) (7 + 1) - m(l - m)( 7 - 2) (7 - 3)] if n > m, 



( 7 +l) 7 ( 7 -l)( 7 -2)( 7 -3) 



(7 + 7 2 + 6772 + 6n 2 ) 



if n = m, 



2( 7 +l) 7 (7-l)(7-2)( 7 -3) V n^lYn 

[(2 + m)(l + m)7(7 + 1) - 2n(l + m)( 7 - 3) (7 + 1) - n(l - n)(7 - 2) (7 - 3)] if n< m. 

(3.15) 

We can derive similar expressions for all even integers beyond 6, i.e. a = 8, 10, . . . where 
we have, for n>m, that 



(-l) n+m p 



(7)n r( 7 ) 



n +m ^ (f)» r(7~f) , / (7)n(7)m ^ (-m).( 7 ~ f ).(! ~ f ), 

n!m! s!( 7 ) a (l - f - n) 



(3.16) 



and the matrix elements with < n < m are incorporated by using the symmetry 
property, i.e. r m " = . For N -dimensional case, the matrix elements of the singular 
operator r~ a can be easily found in analogy with (3.10) 



r -a = (_!)«+™^f (f )" r ^ 7Ar ~ t) . /(7Jv)n(7Jv)r 



n\m\ 



, -m, 7 ^-f,l-f 
7jv,l-n-f 



(3.17) 



(7)n r(7Ar) 

The results for the special cases a = 2,4,6, . .. can be obtained in a similar fashion to 



3.14-3.16 through the substitution of 7 by j N = 1 + + (Z + f - l) 2 . 
3.2 Matrix elements of the power-law potentials r q : 

We now use the orthonormal eigenfunctions (3.3) to compute the matrix elements for 

the power-law potential operators r q , q = 2,4, 6, In analogy with (3.10), this can 

be achieve by means of the identity: for m and n no n- negative integers and + q > 



00 



,2 7+q -l e -0r' iFi 



(3r 2 ) 1 F 1 



- rn 



Pr 2 )dr 



(-|)nr(7 + |) 



, -m, 7 +I,l + I 
3 ^ 21 7,1 + I ~n 



(3.18) 

In the case of q being a positive even number ( q = 2,4, 6, . . . ), the Clausen hypergeo- 
metric function 3^2 in (3.18) can be further simplified. Indeed in this case, we prove 
the following result. 



Theorem 1: For t = 0, 1, 2, . . . , | and q = 2, 4, 6, . . . , the matrix elements of the power- 
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law potential r g , q = 2,4, 6, . .. in terms of the orthonormal functions (3.3) are 



■2 V 



(7+§)™r( 7 +f)(m+§)! 



fit V (l)mT(-y)m\ 

r (7+f )(7+m-<) 9 (m-t+1) q 
2 2_ 

(-i)*-^5t!r(7)( 7 ). 



t 

(7) m+ 3-t(™+*-*)! J? " 
r(7+f) (l-7-m)^ (-m)| I ( _ ^ (7+m)j (m+1) . 

( 7+m -2) j ( m+ l_|) :j Jp 

- t 



(-*)j(7+ra+%-t)j(%+m-t+l)j 
(7+m-t) j (m+l-t)j j! 



if n > m + § 



if n = m + § 



ifn = m+ |- t 



(-/3)3 r( 7 ) |! ( 7 ), , =0 

2 J 

r (7+f )(7+n-t) a (n-t+l) 2 



if m = n 



(-!)'-« /m!r( 7 )( 7 )< 



(7) B+ J_ t (n+i-*)l ^ (| 



(-* )j(7+»+f-t)j(f +n-t+l)j 
(7+«-t) J (n+l-t) 3 j! 



ifm = n+ | — t 



32 V 



(7+f).r( 7 +f)(»+f)! 
(7)„r( 7 )n! 



if to = n + § 



if to > n + § 



(3.19) 



Proof: From (3.18), we have 



m+n o—q/2 



(-l) m+ "/3 



r(7+f) / (7)n(7), 
r(7)(7)nV 7i\m\ 



3F2 



-m,7+ f,l 



7,1 + 



and the problem is now reduced to the simplification of the product 



Using the series representation for 3 F 2 , we can write 



3^2 



-m,7 + |,l + | 
7,1 + 1 -n 



E 



(-m) fc (7 +§)*(-§-*), 



where we have used the identity (a - n) k = (1 - a)„(a) fe /(l — a — k) n . From the definition 
of the Pochhammer symbol (-§ - fe)„ we have that 7 = for n > m + § and for 
n = I + to we have 

J _ (-l) f (7+fUm+§)! 

(7)m 



Finally, for n = m + § - t , i = 0, 1, 2, . . . , § , 
(-m) fc (7+§) fc (-f-A), 



fc=0 



(7)fc 



- E 

k—m — t 



(-m) fc ( 7 +§)*(-§-*) 



m+%—t 



(7)* *;! 



The completion of the proof then follows by shifting the index j = k -m + t of the 
finite sum. □ 
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As consequence of Theorem 1, for q = 2 and 7 > -1 we have that 



ro, 



if n > to + 1 



( 9~ 1 x /(m + 1)(7 + to), ifn = m+l 



/3" 1 (7 + 2n), 



/ 3" 1 v / (" + 1 )(7 + "), ifm = n + l 



if n = m 



(3.20) 



if to > n + 1. 



Furthermore, the explicit formula for the matrix elements of the operator r 4 now 
reads 



if n > to + 2 



/3~ 2 ^/(m + 1)(to + 2) (7 + m)(7 + to + 1), if n = to + 2 

2/3" 2 (7 + 2m + 1)v/(to + 1)(7 + m), if n = to + 1 

/3~ 2 (7 + 6m 2 + 67TO + 7 2 ), if n = to 

2/T 2 ( 7 + 2n + l),/(n + 1)( 7 + n), if to = n + 1 

/3- 2 v /(n + \){n + 2) (7 + n)( 7 + n + 1), if to = n + 2 

I 0, if m > n + 1. 



(3.21) 
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4. Applications and Numerical Results 

4.1 Spiked Harmonic oscillator Hamiltonians 

There are several reasons for the interest in the spiked harmonic oscillator Hamiltonian 
(1.2) and its extension, the so-called generalized spiked harmonic oscillator, 

+ + £ + A. (4,) 

First, it represents the simplest model of certain realistic interaction potentials in 
atomic, molecular and nuclear physics, and second, its interesting intrinsic properties 
from the viewpoint of mathematical physics: (1) an eigenvalue of the perturbed op- 
erator may not converges to the original one as A — > (the Klauder Phenomenon) 
and (2) the perturbation series is ordered in fractional powers of A , and in the cases 
a > 5/2 the regular Rayleigh-Schrodinger perturbation theory fails badly. 

We shall consider the problem initially in TV = 3 spatial dimensions. It was proven 
earlier [44] that the set of L 2 [0,oo) -functions {ip n (r)}%Lo as defined by (3.3), is a 
complete orthonormal basis for the Hilbert space L 2 [0,oo) . This basis was starting 
point for perturbative expansions and variational analysis of the Hamiltonians (1.2) 
and (4.1). The main approach of the earlier variational investigation of (1.2) was the 
re-writing of the Hamiltonian as 

»=-^+n + ^ + (--^y (4.2) 

The parameter A serves as an extra degree of freedom that can be used to accelerate 
the convergence to the exact eigenvalues through the minimization of the eigenvalues 
of the diagonalizable D x D symmetric matrix. Straightforward calculations using 
(3.11) and (3.13) show that the matrix elements H mn of the Hamiltonian (4.2) are 
( m,n = 0, 1, 2, .., D — 1 , n>m) 



H mn -2(2n + 1 )5 nm + (-l) A— ^-^—^3^ 7;1 _^_ f 2 



( 1 yn+n (7 -§)(7- |) / nl(7) ro 

1 ' 7-1 V m! (7)n : 



(4.3) 



where 27 > a and the matrix elements with < n < m are incorporated by using the 
symmetry property of the matrix. In order to apply the method discussed in Section 
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2, we write (1.2) in more extended form 

d 2 „ o A „s9^A A 



H = -i*+ Br2 + ^ + (l-Br+{^-^), (4-4) 

In this case, the matrix elements of the Hamiltonian (4.4) assume the form 

H mn = 2/3(2n + 7 )J ron + (1 - p 2 )r 2 mn + Xr m a n - Ar~ 2 nl (/3 = VB, m, n = 0, 1, 2, . . . , D - 1) (4.5) 

where r 2 mn is given by (3.20), r m " given by (3.11), and r m 2 n is given by (3.13). In order 
to illustrate the difference between using the expressions (4.3) and (4.5), we restrict 
our calculation to a = 4 . The first variational approximation (subspace of dimension 
1) of the ground-state eigenvalues of the spiked harmonic oscillator Hamiltonian yields 
by means of (4.3) the approximation 

e = min {7+I + - ^ 7 + ^^ — A, (l = 1 + -VTTIX > 2) (4.6) 

° A>0.75 \ ' (7-l)( 7 -2) 4(7-1)/' W 2 ' K ' 

while the matrix elements (4.5) yield the approximation 



eo = nun 



"'•13+-, ^ 7z+ At P Tt I> (7=l + ^v / T+4l>2,/3 = v A B). (4.7) 



A>0.75,B>0 \/3 (7-l)(7-2) 4(7-1) J' W 2 

The minimization of (4.6) over the parameter A yields 

e = I + I ^8A — 1 + 4V4A2 - A^J ' + ^6A - - 3^4 A 2 - ^8A - 1 + 4 V / 4A2 _ A^ 7 

which implies the upper bound e = 21.427 793 for A = 1000; while the minimization 
of (4.7) over A and B yields e = 21.374 087 with first decimal place exact. In Table 
I, we present a comparison between eigenvalues computation using (4.6) and (4.7) 
for A = 0.1 to illustrate the increased in the rate of convergence obtained when using 
our new approach. It should be noted that the optimization over the parameter A 
of the 1000 x 1000 -diagonalizable matrix yields an upper bound of E A = 3.575 557 
with A « 6.076 , while the minimization and diagonalization of the symmetric matrix 
over the parameters A and B greatly reduce the number of the matrix elements 
needed by a ratio of approximately 10 : 1 . As shown in Table I, 100 x 100 matrix 
is sufficient to achieve an exact eigenvalues of E A - B = 3.575 552 . Because we have 
established simple formulas for the matrix elements in the cases a = 4 and a = 6 
(given by (3.14) and (3.15) respectively), the determination of the energy values to 
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any desired accuracy reduced to an easy task as indicated in Table II where we 
report our eigenvalue computation for the case a = 6 and for different values of the 
parameter A . A heuristic scheme for ascertaining the eigenvalues to any required 
number of digits is as follows. The eigenvalues obtained from successive levels, such 
as ( 1 x 1, 2 x 2, . . . ), of the truncated matrix are compared, and the calculation ceases 
when the successive eigenvalue agree with each other up to the prescribed decimal 
place. Further advantage of the variational approach presented here is the amount 
of information that we get about the spectrum of the Hamiltonian every time we 
compute the eigenvalues via the diagonalization and minimization. Indeed, we obtain, 
for D x D -matrix, a set of upper bounds for the eigenvalues E ,E 1 , . . . , Eu-\ . Each can 
be improved by either an increase in the dimension of the matrix, or by extracting the 
desired level through the diagonalization and subsequent minimization with respect 
to parameters A and B . For the N -dimensional case, the matrix elements of the 
singular operator r~ a turn out to be 

<lUr l^>-(-l) r(w) V n!m , 3^ lNA _ q __ n 1J. (4.8) 

Matrix elements for the special cases of a = 2, 4, 6, . . . are obtained by substituting 
in Eqs.(3. 13-15) for 7 the expression 7^ , where 7^ = 1 + \J A + (A + ±) 2 . The matrix 
elements of the spiked harmonic oscillator Hamiltonian now turn out to be very similar 
to those in Eq.(4.5), namely 

H mn = 2[3(2n + lN )5 mn + (1 - B)r 2 rnn + \r m a n - Ar m 2 n . (4.9) 

In Table III, upper bounds , obtained by the optimization of the eigenvalues of 
a 10 x 10 -matrix over the parameters A and B . The results are reported for the 
Hamiltonian H = + A( ^ 1) + r 2 + ^ffl where A = i(A + 2Z-3) for dimension N = 2 
to 10 with the angular momentum 1 = 0. 

4.2 Anharmonic Singular Hamiltonian 

The higher-order anharmonic singular Hamiltonians (1.3) have attracted much atten- 
tion recently [47-66]. This is in part because the study of the relevant Schrodinger 
equation with anharmonic potentials provides understanding and insight for the cor- 
responding physical problems, and also because the determination of its energy is 
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itself a challenging problem. In 3-dimensional space, there are two main methods 
for dealing with the anharmonic potentials V(r) = ar 2 + br~ i + cr~ 6 . A method due 
to Varshni [56] is based on an ansatz for the eigenfunctions, sufficient conditions on 
parameters to yield exact solutions, and a limit from initial box confinement. The 
other method, mainly proposed by Znojil [47-48], relies on a Laurent series ansatz for 
the eigenfunctions, which converts the Schrodinger equation into a difference equa- 
tion which is solved by the use of continued fractions. An interesting study related to 
Varshni's idea [56] for the potential V(r) in 2-dimensions was proposed recently by 
Shi-Hai Dong and Zhong-Qi Ma [60]. 

The method discussed in Section 2 of the present article provides a uniformly 
simple, straightforward and very efficient way of yielding accurate energies of the 
entire spectrum of the anharmonic potentials V(r) not only in 1 or 2 dimensions 
but actually in arbitrary dimensions with arbitrary angular momentum number I = 
0, 1,2, We start with a wider class of anharmonic singular Hamiltonian given by 

d 2 A (A + 1) 9 <J 2 &3 &4 , . „. 

H = -^2+^^+ a ^ 2 + ^ + ^ + ^ ( 4 - 10 ) 

where A = |(M - 3) , and M = N + 21 . Clearly, the case V(r) = air 2 + a 3 r~ 4 + a 4 r~ 6 
appears as special case with a 2 — . Following the procedure discussed in Section 2, 
we write the Hamiltonian as 

d 2 , A(A + 1) + A, . 9 (a 2 - A) a 3 

The matrix elements of the Hamiltonian (4.11) take the form ( m,n = 0, 1,2, ... ) 

H mn = 2VS(2n + lN )5 mn + (ai - B)r 2 mn + (a 2 - A)r m 2 n + a 3 r m 4 n + a 4 r m 6 „ (4.12) 

for 7 > 3 . To analyze the precision of the method proposed here, we compare our 
results with some special cases for which the exact eigenvalues are known. The case 
of a\ = a 3 = a 4 = 1, a 2 = 0, which yields the ground-state energy E = 5 has been 
analyzed by Znojil [49], Guardiola and Ros [54], and Buendia et al [51] by different 
techniques. Table (IV) shows the exact eigenvalues of E = 5 can be reached with 
the diagonalization of 50 x 50 -matrix. It should be noted however we have fixed the 
dimension of the matrix to 50 x 50 , but the particular value E = 5 can be reached with 
far fewer matrix elements, indeed a 30x30 matrix is sufficient to achieve such accuracy. 
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Further, the exact energies of 7,7,11 corresponding to (ai,a 3 ,a 4 ) = (1,9,9), (1,-7,49) , 
and (1,45,225) respectively follow simply with the optimization of the diagonalizable 
40 x 40 -matrix with (A, B) = (17.47, 5.69), (18.86, 5.53) , and (17.92,5.40). These results 
simply indicate the generality and the efficient of our approach. Note, the case of 
(01,03,04) = (1,-7,49) also reflect the applicability of the method in the case of the 
parameter a 3 is negative. It is quite clear from Tables (IV) and (V) the generality of 
the method proposed here. In Table (V), we illustrate the applicability of the method 
to the problem of obtaining the energies in different dimensions. Similar results for 
different excited states can be easily reproduce. All the eigenvalues quoted in Tables 
I-V agree with the numerical solutions of the correspondence Schrodinger equation. 
Generally speaking, the precision of the energies to any number of decimal places can 
be easily achieve by increasing the dimension of the matrix. 

5. Conclusion 

We have developed an effective variational method to study a large family of singular 
Hamiltonians. A key feature of this work is the establishment and simplification of 
closed-form analytical expressions for the matrix elements with respect to a basis 
derived from a soluble singular problem. These formulas are general in the sense 
that they include two pameters from the basis which can then be used to optimize 
the matrix eigenvalues obtained for the problem in hand. The improved variational 
approach yields faster energy convergence than was possible earlier. 
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Table (I): Upper bounds E A for H = —A + r 2 + — ^ are obtained by diagonalization then 
minimization of the D x D matrix, only over the parameter A . E A ' B arc the corresponding values 
minimized over both A and B. The eigenvalue 3.575 552 (exact to 7 places) can be easily verified 
by direct numerical integration of Schrodinger's equation. 



D x D 


E A 




1 x 1 


3.745 811 
(A w 1.52) 


3.664 281 
(Aw 1.92,5 w 1.62) 


10 x 10 


3.602 189 
(A « 2.84) 


3.582 194 
(A w 4.75,5 w 11.66) 


20 x 20 


3.588 143 
(A w 3.68) 


3.576 773 
(Aw 8.14,5 w 32.22) 


100 x 100 


3.577 007 
(A w 7.44) 


3.575 552 
(A w 9.73, 5 w 297.23) 


200 x 200 


3.576 015 
(A w 10.39) 


3.575 552 
(A w 3.76, 5 w 873.58) 
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Table (II): Upper bounds E A for H = —A + r 2 + obtained by diagonalization followed 
by minimization of the D x D matrix, only over the parameter A . E A ' B are the corresponding 
values minimized over both A and B. The results are displayed for different values of the potential 
coupling A . The eigenvalues are correct for the 7 digits, as can be easily verify by direct numerical 
integration of Schrodingcr's equation. 



A 


E A 


E A,B 


1000 


12.718 617 


12.718 617 






(32 x 32, Am 20.52) 


(15 x 15, A m 54.41,5 


m 4.51) 


100 


8.413 358 


8.413 358 






(65 x 65, A m 19.61) 


(22 x 22, Am 8.88, B t 


a 9.76) 


10 


6.003 209 


6.003 209 






(150 x 150, A m 9.71) 


(30 x 30, A m 9.66, B » 


i 21.41) 


1 


4.659 940 


4.659 940 






(350 x 350, A m 17.79) 


(45 x 45, Am 18.34, B t 


a 56.77) 


0.1 


3.915 665 


3.915 665 






(1000 x 1000) 


(80 x 80, A m 4.47, B m 


176.63) 


0.01 


3.505 492 


3.505 455 






(1000 x 1000) 


(100 x 100, Am 20.88,5 


m 348.92) 
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Table (III) Upper bounds for H = + ^±1) + r 2 + iroo for dimcnsion N = 2 to 
10 , obtained by diagonalization then minimization of the 10 x 10 matrix over A and B . 



N 

1 V 




2 


21 350 246 (A ?a 71 44 B«l 775") 


3 


21.369 463 (A w 71.27, B w 1.774) 


4 


21.427 056 (A w 70.79, S w 1.772) 


5 


21.522 860 (A w 69.96, B w 1.769) 


6 


21.656 596 (A w 68.82, B w 1.764) 


7 


21.827 883 (A w 67.34, S w 1.757) 


8 


22.036 232 (A w 65.45, B w 1.749) 


9 


22.281 057 (A w 63.35, S w 1.740) 


10 


22.561 680 (A w 60.81, S w 1.726) 
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Table (IV): Upper bounds for the Hamiltonian V(r) = ar 2 + br~ 4 + cr~ 6 for different values 
of the parameters a, b, and c . The eigenvalues are exact for the 7 digits shown, as confirmed 
numerically (or, for the first row, known exactly). 



a 


b 


c 


E u 


1 


1 


1 


5.000 000 
(50 x 50, Aw 25.51, B « 54.35) 


1 


10 


1 


6.679 054 
(50 x 50, Ak 9.65, B w 32.24) 


1 


1 


10 


6.140 123 
(50 x 50, Ak 5.56, B w 38.67) 


1 


10 


10 


7.138 261 
(50 x 50, A w 7.83, B w 31.94) 


1 


too 


100 


11.791 771 

(50 x 50, A w 5.15, B w 5.00) 


1 


1000 


1000 


21.885 192 
(50 x 50, A w 19.42, B w 5.00) 
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Table (V): Upper bounds E N for H = + + r 2 + £ + ^ for dimension N = 2 

to 10 , obtained by diagonalization then minimization with respect A and B of a 30 x 30 matrix. 



N 

1 V 




2 


1 2 704 404 (A ss 5 05 B w 5 53") 


3 


12.735 264 (A w 11.73, B w 6.40) 


4 


12.827 666 (A w 5.05, S w 7.82) 


5 


12.981 081 (A w 5.97, B w 7.26) 


6 


13.194 635 (A w 6.25, S w 7.28) 


7 


13.467 115 (A w 6.12, B w 7.35) 


8 


13.796 990 (A w 6.16, B w 7.46) 


9 


14.182 423 (A w 6.31, S w 7.32) 


10 


14.621 300 (A w 3.11, S w 7.23) 
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